coord<-read.csv("coord_data.csv")

### TO CREATE ADJACENCY MATRIX ##
cr <- paste(coord$col, coord$row)

####### ADJACENCY MATRIX 1 #######
###ADJ1C1
adj1c1.cr<- paste(coord$col-1, coord$row-1)
coord$adj1c1<- coord$gid[match(adj1c1.cr,cr)]

###ADJ1C2
adj1c2.cr<- paste(coord$col, coord$row-1)
coord$adj1c2<- coord$gid[match(adj1c2.cr,cr)]

###ADJ1C3
adj1c3.cr<- paste(coord$col+1, coord$row-1)
coord$adj1c3<- coord$gid[match(adj1c3.cr,cr)]

###ADJ1C4
adj1c4.cr<- paste(coord$col-1, coord$row)
coord$adj1c4<- coord$gid[match(adj1c4.cr,cr)]

###ADJ1C5
adj1c5.cr<- paste(coord$col+1, coord$row)
coord$adj1c5<- coord$gid[match(adj1c5.cr,cr)]

###ADJ1C6
adj1c6.cr<- paste(coord$col-1, coord$row+1)
coord$adj1c6<- coord$gid[match(adj1c6.cr,cr)]

###ADJ1C7
adj1c7.cr<- paste(coord$col,coord$row+1)
coord$adj1c7<- coord$gid[match(adj1c7.cr,cr)]

###ADJ1C8
adj1c8.cr<- paste(coord$col+1,coord$row+1)
coord$adj1c8<- coord$gid[match(adj1c8.cr,cr)]

####### ADJACENCY MATRIX 2 #######
###ADJ2C1
adj2c1.cr<- paste(coord$col-2, coord$row-2)
coord$adj2c1<- coord$gid[match(adj2c1.cr,cr)]

###ADJ2C2
adj2c2.cr<- paste(coord$col-1, coord$row-2)
coord$adj2c2<- coord$gid[match(adj2c2.cr,cr)]

###ADJ2C3
adj2c3.cr<- paste(coord$col, coord$row-2)
coord$adj2c3<- coord$gid[match(adj2c3.cr,cr)]

###ADJ2C4
adj2c4.cr<- paste(coord$col+1, coord$row-2)
coord$adj2c4<- coord$gid[match(adj2c4.cr,cr)]

###ADJ2C5
adj2c5.cr<- paste(coord$col+2, coord$row-2)
coord$adj2c5<- coord$gid[match(adj2c5.cr,cr)]

###ADJ2C6
adj2c6.cr<- paste(coord$col-2, coord$row-1)
coord$adj2c6<- coord$gid[match(adj2c6.cr,cr)]

###ADJ2C7
adj2c7.cr<- paste(coord$col+2, coord$row-1)
coord$adj2c7<- coord$gid[match(adj2c7.cr,cr)]

###ADJ2C8
adj2c8.cr<- paste(coord$col-2, coord$row)
coord$adj2c8<- coord$gid[match(adj2c8.cr,cr)]

###ADJ2C9
adj2c9.cr<- paste(coord$col+2, coord$row)
coord$adj2c9<- coord$gid[match(adj2c9.cr,cr)]

###ADJ2C10
adj2c10.cr<- paste(coord$col-2, coord$row+1)
coord$adj2c10<- coord$gid[match(adj2c10.cr,cr)]

###ADJ2C11
adj2c11.cr<- paste(coord$col+2, coord$row+1)
coord$adj2c11<- coord$gid[match(adj2c11.cr,cr)]

###ADJ2C12
adj2c12.cr<- paste(coord$col-2, coord$row+2)
coord$adj2c12<- coord$gid[match(adj2c12.cr,cr)]

###ADJ2C13
adj2c13.cr<- paste(coord$col-1, coord$row+2)
coord$adj2c13<- coord$gid[match(adj2c13.cr,cr)]

###ADJ2C14
adj2c14.cr<- paste(coord$col, coord$row+2)
coord$adj2c14<- coord$gid[match(adj2c14.cr,cr)]

###ADJ2C15
adj2c15.cr<- paste(coord$col+1, coord$row+2)
coord$adj2c15<- coord$gid[match(adj2c15.cr,cr)]

###ADJ2C16
adj2c16.cr<- paste(coord$col+2, coord$row+2)
coord$adj2c16<- coord$gid[match(adj2c16.cr,cr)]


### Remove excess vectors from the workspace ######
rm(list = setdiff(ls(), "coord"))


### TO EXTRACT NEIGHBORS' CONFLICT VALUE ##

### MERGE DATA
rdata <- read.delim("Replication Data.tab", header=T)
load("cell_data.RData")
data_all <- merge(rdata, cell_data, by = "gid")
data<-merge(data_all, coord, by=c("gid", "cell_dum_07", "cmean_07"))

rm(rdata)
rm(cell_data)
rm(data_all)

###V_ADJ1C1 through V_ADJ2C16###
data$V_adj1c1<- data$conf2008_dum[match(data$gid,data$adj1c1)]
data$V_adj1c2<- data$conf2008_dum[match(data$gid,data$adj1c2)]
data$V_adj1c3<- data$conf2008_dum[match(data$gid,data$adj1c3)]
data$V_adj1c4<- data$conf2008_dum[match(data$gid,data$adj1c4)]
data$V_adj1c5<- data$conf2008_dum[match(data$gid,data$adj1c5)]
data$V_adj1c6<- data$conf2008_dum[match(data$gid,data$adj1c6)]
data$V_adj1c7<- data$conf2008_dum[match(data$gid,data$adj1c7)]
data$V_adj1c8<- data$conf2008_dum[match(data$gid,data$adj1c8)]

data$V_adj2c1<- data$conf2008_dum[match(data$gid,data$adj2c1)]
data$V_adj2c2<- data$conf2008_dum[match(data$gid,data$adj2c2)]
data$V_adj2c3<- data$conf2008_dum[match(data$gid,data$adj2c3)]
data$V_adj2c4<- data$conf2008_dum[match(data$gid,data$adj2c4)]
data$V_adj2c5<- data$conf2008_dum[match(data$gid,data$adj2c5)]
data$V_adj2c6<- data$conf2008_dum[match(data$gid,data$adj2c6)]
data$V_adj2c7<- data$conf2008_dum[match(data$gid,data$adj2c7)]
data$V_adj2c8<- data$conf2008_dum[match(data$gid,data$adj2c8)]
data$V_adj2c9<- data$conf2008_dum[match(data$gid,data$adj2c9)]
data$V_adj2c10<- data$conf2008_dum[match(data$gid,data$adj2c10)]
data$V_adj2c11<- data$conf2008_dum[match(data$gid,data$adj2c11)]
data$V_adj2c12<- data$conf2008_dum[match(data$gid,data$adj2c12)]
data$V_adj2c13<- data$conf2008_dum[match(data$gid,data$adj2c13)]
data$V_adj2c14<- data$conf2008_dum[match(data$gid,data$adj2c14)]
data$V_adj2c15<- data$conf2008_dum[match(data$gid,data$adj2c15)]
data$V_adj2c16<- data$conf2008_dum[match(data$gid,data$adj2c16)]

### TO EXTRACT NEIGHBORS' CELLPHONE VALUE ##
data$C_adj1c1<- data$cell_dum_07[match(data$gid,data$adj1c1)]
data$C_adj1c2<- data$cell_dum_07[match(data$gid,data$adj1c2)]
data$C_adj1c3<- data$cell_dum_07[match(data$gid,data$adj1c3)]
data$C_adj1c4<- data$cell_dum_07[match(data$gid,data$adj1c4)]
data$C_adj1c5<- data$cell_dum_07[match(data$gid,data$adj1c5)]
data$C_adj1c6<- data$cell_dum_07[match(data$gid,data$adj1c6)]
data$C_adj1c7<- data$cell_dum_07[match(data$gid,data$adj1c7)]
data$C_adj1c8<- data$cell_dum_07[match(data$gid,data$adj1c8)]

data$C_adj2c1<- data$cell_dum_07[match(data$gid,data$adj2c1)]
data$C_adj2c2<- data$cell_dum_07[match(data$gid,data$adj2c2)]
data$C_adj2c3<- data$cell_dum_07[match(data$gid,data$adj2c3)]
data$C_adj2c4<- data$cell_dum_07[match(data$gid,data$adj2c4)]
data$C_adj2c5<- data$cell_dum_07[match(data$gid,data$adj2c5)]
data$C_adj2c6<- data$cell_dum_07[match(data$gid,data$adj2c6)]
data$C_adj2c7<- data$cell_dum_07[match(data$gid,data$adj2c7)]
data$C_adj2c8<- data$cell_dum_07[match(data$gid,data$adj2c8)]
data$C_adj2c9<- data$cell_dum_07[match(data$gid,data$adj2c9)]
data$C_adj2c10<- data$cell_dum_07[match(data$gid,data$adj2c10)]
data$C_adj2c11<- data$cell_dum_07[match(data$gid,data$adj2c11)]
data$C_adj2c12<- data$cell_dum_07[match(data$gid,data$adj2c12)]
data$C_adj2c13<- data$cell_dum_07[match(data$gid,data$adj2c13)]
data$C_adj2c14<- data$cell_dum_07[match(data$gid,data$adj2c14)]
data$C_adj2c15<- data$cell_dum_07[match(data$gid,data$adj2c15)]
data$C_adj2c16<- data$cell_dum_07[match(data$gid,data$adj2c16)]

### SUM OF VIOLENCE ###
violence1<-data[,c(66:73)]
violence2<-data[,c(74:89)]


data$V_adj1sum<-rowSums(violence1, na.rm=TRUE)

data$V_adj1bin <- 0
for (i in 1:10674){ 
  if (data$V_adj1sum[i]>=1) {
    data$V_adj1bin[i]<-1}
}

data$V_adj2sum<-rowSums(violence2, na.rm=TRUE)

data$V_adj2bin <- 0
for (i in 1:10674){ 
  if (data$V_adj2sum[i]>=1) {
    data$V_adj2bin[i]<-1}
}

### MATRIX MULTIPLICATION: VIOLENCE *T(CELL) ###
cell1<-data[,c(90:97)]
cell2<-data[,c(98:113)]
data$inter_VC1<-rowSums(cell1*violence1, na.rm=TRUE)
data$inter_VC2<-rowSums(cell2*violence2, na.rm=TRUE)

## Triple interaction: neighbor violence and cell with ref cell phone #####
data$triple_VC1 <- data$inter_VC1 * data$cell_dum_07
data$triple_VC2 <- data$inter_VC2 * data$cell_dum_07

###########################################################################

## Models, tables from the paper #####

## Requires the older version of Zelig for the robust 
## standard errors to work

# Make sure that Zelig_3.5.5.tar.gz is in your working directory
install.packages("Zelig_3.5.5.tar.gz", repos=NULL, type="source")
require(Zelig)

## Table 1 ######################

# Logit model from P & K, no ref cell phones ########
m1_ph <- zelig(conf2008_dum ~ pre2000_count + bdist1 + capdist + 
                    pop2005 + mnt + irri + gcppc00,
                  model="logit",
                  data = data,robust=TRUE)
summary(m1)

# Logit model from P&K, with ref cell phones #########
m2_ph <- zelig(conf2008_dum ~ pre2000_count + bdist1 + capdist + 
              pop2005 + mnt + irri + gcppc00 + cell_dum_07,
            model="logit",
            data = data,robust=TRUE)
summary(m2)

# Re-logit model from P&K, with ref cell phones #########
m3_ph <- zelig(conf2008_dum ~ pre2000_count + bdist1 + capdist + 
              pop2005 + mnt + irri + gcppc00 + cell_dum_07,
            model="relogit",
            data = data,robust=TRUE)
summary(m4)

# Fixed-effects model from P&K, with ref cell phones ######
m4_ph <- zelig(conf2008_dum ~ pre2000_count + bdist1 + capdist + 
              pop2005 + mnt + irri + gcppc00 + cell_dum_07 +
              factor(cow),
            model="ls", data=data, robust=TRUE)
summary(m4)

# Table for P&K models ############
stargazer(m1_ph, m2_ph, m3_ph, m4_ph,
          keep = c("Intercept", "pre2000_count", "bdist1", "capdist", "pop2005", 
                   "mnt", "irri", "gcppc00", "cell_dum_07"))



###############################################################

## Tables for Spatial Models in Section 3 ######

# Logit model, cell phones, first-neighbor violence #####
m2_first <- zelig(conf2008_dum~pre2000_count+bdist1+capdist+
                    pop2005+mnt+irri+gcppc00+cell_dum_07 + 
                    V_adj1sum,
                  model="logit",
                  data=data,robust=TRUE)
summary(m2_first)

# Logit model, cell phones, second-neighbor violence #####
m2_second <- zelig(conf2008_dum~pre2000_count+bdist1+capdist+
                     pop2005+mnt+irri+gcppc00+cell_dum_07 + 
                     V_adj2sum,
                   model="logit",
                   data=data,robust=TRUE)
summary(m2_second)

# Logit model, cell phones, both neighbor violence ######
m2_both <- zelig(conf2008_dum~pre2000_count+bdist1+capdist+
                     pop2005+mnt+irri+gcppc00+cell_dum_07 + 
                     V_adj1sum + V_adj2sum,
                   model="logit",
                   data=data,robust=TRUE)
summary(m2_both)

# Table of neighbor violence (Table 2) #######
stargazer(m2_first, m2_second, m2_both,
          keep = c("Intercept", "pre2000_count", "bdist1", "capdist", "pop2005", 
                   "mnt", "irri", "gcppc00", "cell_dum_07", "V_adj1sum", "V_adj2sum"))

### FIRST DIFFERENCES AND PREDICTED VALUES ###
library(mvtnorm)

coef <- coef(m2_first) 
varcov <- vcov(m2_first) 


neighbors1_coef<-cbind(pre2000_count, bdist1, capdist, pop2005, 
                       mnt, irri, gcppc00, cell_dum_07,V_adj1sum)
neighbors1_vec<- apply(X = cbind(1,neighbors1_coef), MARGIN = 2, FUN = mean, na.rm=T) 

set.seed(1234)
sim.betas <- rmvnorm(n = 10000,
                     mean = coef,
                     sigma = varcov)

presvals <- 0:8
ev.ests <- matrix(data = NA, ncol = length(presvals),
                  nrow=10000)

for(j in 1:length(presvals)){
  vec.new <- neighbors1_vec
  vec.new[10] <- presvals[j]
  ev.ests[,j] <- pnorm(vec.new%*%t(sim.betas))
}
plot(presvals, apply(ev.ests,2,mean), ylim=c(0,1),
     xlab = "Number of neighbor cells with violence",
     ylab = "Probability of violence in the reference cell")
segments(x0 = presvals, x1 = presvals,
         y0 = apply(ev.ests, 2, quantile, .025),
         y1 = apply(ev.ests, 2, quantile, .975)
)

################################################################

## The triple interaction models ########

# Logit model, first neighbor, with triple interaction #########
m1_first_triple <- zelig(conf2008_dum~pre2000_count+bdist1+capdist+
                           pop2005+mnt+irri+gcppc00 + cell_dum_07 +
                           V_adj1bin + triple_VC1
                         ,model="logit",
                         data=data,robust=TRUE)
summary(m1_first_triple)

# Logit model, second neighbor, with triple interaction #######
m1_second_triple <- zelig(conf2008_dum~pre2000_count+bdist1+capdist+
                            pop2005+mnt+irri+gcppc00 + cell_dum_07 +
                            V_adj2bin + triple_VC2
                          ,model="logit",
                          data=data,robust=TRUE)
summary(m1_second_triple)

# Logit model, both, with triple interaction #######
m1_both_triple <- zelig(conf2008_dum~pre2000_count+bdist1+capdist+
                            pop2005+mnt+irri+gcppc00 + cell_dum_07 +
                            V_adj1bin + V_adj2bin + triple_VC2
                          ,model="logit",
                          data=data,robust=TRUE)
summary(m1_both_triple)

## Table 3: Logit models with triple interaction #############
stargazer(m1_first_triple, m1_second_triple, m1_both_triple)
